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Abstract 


Space— Time  Interface- Tracking  Computations  with  Contact 

Between  Solid  Surfaces 

by 

Austin  J.  Buscher 


To  address  the  computational  challenges  associated  with  contact  between  mov¬ 
ing  solid  surfaces,  such  as  those  in  cardiovascular  fluid-structure  interaction  (FSI), 
parachute  FSI,  and  flapping-wing  aerodynamics,  we  introduce  a  space-time  (ST) 
interface-tracking  method  that  can  deal  with  topology  change  (TC).  In  cardiovas¬ 
cular  FSI,  our  primary  target  is  heart  valves.  The  method  is  a  new  version  of  the 
Deforming-Spatial-Domain/Stabilized  ST  (DSD/SST)  method,  and  we  call  it  ST-TC. 
It  includes  a  master-slave  system  that  maintains  the  connectivity  of  the  “parent” 
mesh  when  there  is  contact  between  the  moving  interfaces.  It  is  an  efficient,  practical 
alternative  to  using  unstructured  ST  meshes,  but  without  giving  up  on  the  accurate 
representation  of  the  interface  or  consistent  representation  of  the  interface  motion. 
We  explain  the  method  with  conceptual  examples  and  present  2D  and  3D  test  com¬ 
putations  with  models  representative  of  the  classes  of  problems  we  are  targeting. 
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Chapter  1 


Introduction 


The  material  in  this  chapter  is  from  [59].  Flow  problems  with  moving  boundaries 
and  interfaces  (MBI)  are  encountered  frequently  in  engineering  analysis  and  design. 
They  include  fluid-structure  interaction  (FSI),  fluid-object  interaction  (FOI),  fluid- 
particle  interaction  (FPI),  free-surface  and  multi-fluid  flows,  and  flows  with  solid 
surfaces  in  fast,  linear  or  rotational  relative  motion.  These  problems  pose  some  of  the 
most  formidable  computational  challenges.  The  challenges  include  contact  between 
moving  solid  surfaces  and  other  cases  of  topology  change  (TC),  such  as  those  in 
cardiovascular  FSI,  parachute  FSI,  and  flapping-wing  aerodynamics.  A  method  for 
flows  with  MBI  can  be  viewed  as  an  interface-tracking  (moving-mesh)  technique  or 
an  interface-capturing  (nonmoving-mesh)  technique,  or  possibly  a  combination  of  the 
two. 

In  interface-tracking  methods,  as  the  interface  moves  and  the  spatial  domain  oc¬ 
cupied  by  the  fluid  changes  its  shape,  the  mesh  moves  to  accommodate  this  shape 
change  and  to  follow  (i.e.  “track”)  the  interface.  The  Arbitrary  Lagrangian-Eulerian 
(ALE)  finite  element  formulation  [26]  is  the  most  widely  used  moving-mesh  technique, 
with  increased  emphasis  on  FSI  in  recent  years  (see,  for  example,  [41,  78,  7,  30,  33, 
17,  6,  20,  8,  10,  18,  14,  13,  9,  11,  4,  24,  40,  12,  3,  22,  25,  5,  37,  38,  79,  81,  32,  16,  15, 
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31,  82,  29,  23]).  The  Deforming-Spatial-Domain/Stabilized  Space-Time  (DSD/SST) 
method  [65,  68,  69,  67,  70,  54,  56,  16]  is  also  a  general-purpose  moving-mesh  tech¬ 
nique. 

Moving-mesh  methods  require  mesh  update  methods.  Mesh  update  typically  con¬ 
sists  of  moving  the  mesh  for  as  long  as  possible  and  remeshing  as  needed.  With  the  key 
objectives  being  to  maintain  the  element  quality  near  solid  surfaces  and  to  minimize 
frequency  of  remeshing,  a  number  of  advanced  mesh  update  methods  [63,  27,  66,  70] 
were  developed  to  be  used  with  the  DSD/SST  method,  including  those  that  minimize 
the  deformation  of  the  layers  of  small  elements  placed  near  solid  surfaces. 

Over  the  past  20  years  the  DSD/SST  method  has  been  applied  to  some  of  the  most 
challenging  flow  problems  with  MBI.  The  classes  of  problems  solved  include  the  free- 
surface  and  multi-fluid  flows  [65,  69,  63,  62,  66],  FOI  [65,  68,  69,  62],  aerodynamics  of 
flapping  wings  [39,  48,  50],  flows  with  solid  surfaces  in  fast,  linear  or  rotational  relative 
motion  [62,  66,  49,  47,  12,  61],  compressible  flows  [62],  shallow-water  flows  [66,  44], 
FPI  [62,  66],  and  FSI  [39,  28,  64,  71,  70,  34,  72,  73,  43,  76,  74,  52,  35,  77,  75,  54,  53, 
36,  55,  45,  56,  46,  51,  57]. 

As  mentioned  in  [72],  moving  the  fluid  mechanics  mesh  to  track  a  fluid-solid 
interface  enables  us,  at  least  for  interfaces  with  reasonable  geometric  complexity,  to 
control  the  mesh  resolution  near  that  interface  and  obtain  accurate  solutions  in  such 
critical  flow  regions.  As  also  mentioned  in  [72],  sometimes  the  geometric  complexity  of 
the  interface  may  require  a  fluid  mechanics  mesh  that  is  not  affordable  or  not  desirable 
or  just  not  manageable  in  mesh  moving,  and  this  is  one  of  the  most  common  reasons 
given  for  favoring  an  interface-capturing  method.  This  approach  can  be  seen  as  a 
special  case  of  interface  representation  techniques  where  the  interface  geometry  is 
somehow  represented  over  a  nonmoving  fluid  mechanics  mesh,  the  main  point  being 
that  the  fluid  mechanics  mesh  does  not  move  to  track  the  interfaces.  However,  as 
pointed  out  in  [67],  a  consequence  of  the  mesh  not  moving  to  track  the  interface  is 
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that  for  fluid-solid  interfaces,  independent  of  how  accurately  the  interface  geometry 
is  represented,  the  resolution  of  the  boundary  layer  will  be  limited  by  the  resolution 
of  the  fluid  mechanics  mesh  where  the  interface  is. 

For  interfaces  with  reasonable  geometric  complexity,  if  a  moving-mesh  method 
can  be  used  with  a  reasonable  frequency  of  remeshing  (see  [70]  for  various  remeshing 
options),  its  fluid  mechanics  accuracy  near  the  interface  will  be  superior  to  that  of 
an  nonmoving- mesh  method.  As  pointed  out  in  [72],  “while  it  is  understandable  that 
fixed-mesh  methods  become  more  favored  when  the  interface  geometric  complexity 
appears  to  be  too  high  for  a  moving-mesh  method,  we  need  to  remember  that  there 
is  a  difference  between  making  the  problem  computable  and  obtaining  good  fluid 
mechanics  accuracy  near  the  interface.”  Therefore,  as  also  pointed  out  in  [72],  “it 
is  not  difficult  to  imagine  that  if  we  lower  our  expectations  of  good  fluid  mechanics 
accuracy  near  the  interfaces  with  high  geometric  complexity,  we  can  find  a  number 
of  ways  to  make  the  problem  computable  also  with  moving-mesh  methods,  and  can 
still  expect  to  obtain  good  accuracy  near  the  interfaces  with  reasonable  geometric 
complexity.”  Examples  of  that  were  given  in  [72], 

A  robust  moving-mesh  method  with  effective  mesh  update  can  handle  FSI  or 
other  MBI  problems  even  when  the  solid  surfaces  undergo  large  displacements  (see, 
for  example,  FPI  [62,  66]  with  the  number  of  particles  reaching  1,000  [66],  parachute 
FSI  [70,  72,  73,  75,  53,  55,  46,  51,  57],  flapping- wing  aerodynamics  [48,  50],  and 
wind-turbine  rotor  and  tower  aerodynamics  [61].  It  can  handle  FSI  or  other  MBI 
problems  also  even  when  the  solid  surfaces  are  in  near  contact  or  create  near  TC,  if  the 
“nearness”  is  sufficiently  “near”  for  the  purpose  of  solving  the  problem.  Examples  of 
such  problems  are  FPI  with  collision  between  the  particles  [62,  66],  parachute-cluster 
FSI  with  contact  between  the  parachutes  of  the  cluster  [53,  55,  46,  57],  flapping- wing 
aerodynamics  with  the  forewings  and  hindwings  crossing  each  other  very  close  [48,  50], 
and  wind-turbine  rotor  and  tower  aerodynamics  with  the  blades  passing  the  tower 
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close  [61]. 

As  mentioned  in  [16],  one  of  course  recognizes  that  certain  classes  of  interfaces 
(such  as  free-surface  and  two-fluid  flows  with  splashing)  might  be  too  complex  to 
deal  with  an  interface-tracking  technique  and  therefore,  for  all  practical  purposes, 
require  an  interface-capturing  technique.  The  Mixed  Interface-1 Tracking/ Interface- 
Capturing  Technique  (MITICT)  [66]  was  introduced  in  2001  for  computation  of  flow 
problems  that  involve  both  fluid-solid  interfaces  that  can  be  accurately  tracked  with 
a  moving-mesh  method  and  fluid-fluid  interfaces  that  are  too  complex  or  unsteady 
to  be  tracked.  Those  fluid-fluid  interfaces  are  captured  over  the  mesh  tracking  the 
fluid-solid  interfaces.  The  MITICT  was  successfully  tested  in  2D  computations  with 
solid  circles  and  free  surfaces  [2,  19]  and  in  3D  computation  of  ship  hydrodynamics  [3]. 

In  some  MBI  problems  with  contact  between  the  solid  surfaces,  the  “nearness” 
that  can  be  modeled  with  a  moving-mesh  method  without  actually  bringing  the  sur¬ 
faces  into  contact  might  not  be  “near”  enough  for  the  purpose  of  solving  the  problem. 
Cardiovascular  FSI  with  heart  valves,  where  the  flow  has  to  be  completely  blocked 
at  contact,  is  an  example.  The  Fluid-Solid  Interface- Tracking/Interface-Capturing 
Technique  (FSITICT)  [75]  was  motivated  by  such  cardiovascular  FSI  problems.  In 
the  FSITICT,  we  track  the  interface  we  can  with  a  moving  mesh,  and  capture  over 
that  moving  mesh  the  interfaces  we  cannot  track,  specifically  the  interfaces  where  we 
need  to  have  an  actual  contact  between  the  solid  surfaces.  A  specific  application  of 
the  FSITICT  was  presented  in  [80] ,  where  the  ALE  method  is  used  for  interface  track¬ 
ing,  and  a  fully  Eulerian  approach  for  interface  capturing,  with  some  2D  benchmark 
problems  as  test  computations.  This  specific  application  was  extended  in  [80]  to  2D 
FSI  models  with  flapping  and  contact,  where  the  fully  Eulerian  interface-capturing  is 
complemented  with  mesh  adaptivity. 

There  are  many  types  of  nonmoving-mesh  methods  that  can  compute  MBI  prob¬ 
lems  involving  an  actual  contact  between  solid  surfaces  or  other  cases  of  TC.  The 
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immersed-boundary  methods,  X-FEM,  and  particle  methods  are  the  typical  exam¬ 
ples.  Some  of  those  methods  give  up  on  the  accurate  representation  of  the  interface, 
and  most  give  up  on  the  consistent  representation  of  the  interface  motion.  The 
DSD/SST  formulation  does  not  need  to  give  up  on  either,  even  where  we  have  an 
actual  contact  or  some  other  TC,  provided  that  we  can  update  the  mesh  even  there. 
Using  an  ST  mesh  that  is  unstructured  both  in  space  and  time,  as  proposed  for  con¬ 
tact  problems  in  [66],  would  give  us  such  a  mesh  update  option.  However,  that  would 
require  a  fully  unstructured  4D  mesh  generation,  and  that  is  not  easy  in  computing 
real-world  problems. 

We  want  to  address  the  computational  challenges  associated  with  contact  between 
moving  solid  surfaces  and  other  cases  of  TC,  including  those  in  cardiovascular  FSI, 
parachute  FSI,  and  flapping-wing  aerodynamics,  with  the  primary  target  in  cardio¬ 
vascular  FSI  being  heart  valves.  For  this  purpose,  we  introduce  in  this  paper  an  ST 
interface-tracking  method  that  can  deal  with  TC.  It  is  a  new  version  of  the  DSD/SST 
method,  and  we  call  it  ST-TC.  It  is  a  practical  alternative  to  using  unstructured  ST 
meshes,  but  without  giving  up  on  the  accurate  representation  of  the  interface  or  the 
consistent  representation  of  the  interface  motion,  even  where  there  is  an  actual  con¬ 
tact  between  solid  surfaces  or  other  TC.  The  ST-TC  method  is  based  on  special  mesh 
generation  and  update,  and  a  master-slave  system  that  maintains  the  connectivity  of 
the  “parent”  mesh  when  there  is  a  TC. 

In  Chapter  2,  we  provide,  with  two  hypothetical  cases,  a  context  for  TC  and 
explain  the  master-slave  system  and  its  design.  In  Chapter  3,  we  provide  two  con¬ 
ceptual  examples  that  help  us  explain  the  mesh  update  process.  Three  numerical 
examples  are  presented,  a  2D  flapping  pair  in  Chapter  4,  a  3D  Micro  air  vehicle 
(MAV)  in  Chapter  5,  and  an  aortic  valve  model  with  coronary  arteries  in  Chapter  6. 
The  concluding  remarks  are  in  Chapter  7. 
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Chapter  2 


Topology  change 


The  material  in  this  chapter  is  from  [59].  We  consider  two  hypothetical  cases  of  two 
bars  to  provide  a  context  for  TC.  In  the  first  case,  shown  in  Figure  2.1,  the  bars 


Figure  2.1:  Hypothetical  case  of  two  bars  that  are  initially  coinciding,  with  one  hole 
in  the  fluid  mechanics  domain  (top).  Then  the  red  bar  starts  moving  upward,  creating 
a  second  hole  in  the  domain  (bottom). 


are  initially  coinciding,  with  just  one  hole  in  the  fluid  mechanics  domain.  Then  the 
red  bar  starts  moving  upward,  creating  a  second  hole  in  the  domain.  In  the  second 
case,  shown  in  Figure  2.2,  the  bars  are  initially  aligned  with  connected  ends,  again 
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Figure  2.2:  Hypothetical  case  of  two  bars  that  are  are  initially  aligned  with  connected 
ends,  with  one  hole  in  the  domain  (top).  Then  the  red  bar  starts  a  flapping  motion, 
up  (middle)  and  down  (bottom),  creating  a  second  hole  in  the  domain,  except  when 
their  ends  become  connected  periodically  during  the  flapping  motion. 


with  a  single  hole  in  the  domain.  Then  the  red  bar  starts  a  flapping  motion,  up  and 
down,  creating  a  second  hole  in  the  domain,  except  when  their  ends  become  connected 
periodically  during  the  flapping  motion.  When  the  red  bar  is  in  the  upper  position, 
the  part  of  the  domain  below  it  is  connected  to  the  part  of  the  domain  above  the  blue 
bar.  When  the  red  bar  is  in  the  lower  position,  the  part  of  the  domain  above  it  is 
connected  to  the  part  of  the  domain  below  the  blue  bar. 

These  two  cases  are  representatives  of  the  typical  TC  challenges  we  expect  to  see 
in  the  classes  of  MBI  problems  we  are  targeting.  Especially  the  first  case  is  really  not 
possible  to  treat  in  a  consistent  way  without  using  an  ST  method. 


2.1  Master— slave  system 


We  propose  a  very  simple  technique  in  the  ST  context.  Having  a  constraint  between 
nodes  in  a  finite  element  formulation  is  quite  common.  These  constraints  reduce  the 
number  of  unknowns,  but  in  our  implementation  we  delay  that  unknown  elimination 
until  the  iterative  solution  of  the  linear  systems  encountered  at  nonlinear  iterations 
of  a  time  step.  The  iterative  solution  of  the  linear  systems  is  performed  with  reduced 
number  of  unknowns.  The  technique  is  easy  to  manage  in  a  parallel-computing  en¬ 
vironment,  especially  if  the  preconditioner  is  simple  enough.  Typically  we  assign  a 
master  node  to  each  slave  node,  and  we  use  only  the  unknowns  of  the  master  nodes 
in  iterative  solution  of  the  linear  systems. 

We  can  use  different  master-slave  relationships  at  different  time  levels.  This  is  a 
practical  alternative  to,  but  less  general  than,  using  ST  meshes  that  are  unstructured 
in  time.  Still,  we  can  use  this  concept  to  deal  with  the  TC  cases  considered  above,  and 
the  important  point  is  that  the  connectivity  of  the  “parent”  mesh  does  not  change. 
Consequently,  the  distribution  model  in  the  parallel-computing  environment  does  not 
change  during  the  computations  with  moving  meshes. 

With  this  technique,  we  need  to  implement  one  more  functionality.  We  exclude 
certain  elements  from  the  integration  of  the  finite  element  formulation.  The  exclusion 
principles  are  given  below. 

•  Exclude  all  spatial  elements  with  zero  volume  from  the  spatial  integration. 

•  Exclude  all  ST  elements  with  zero  ST  volume  from  the  ST  integration. 

•  We  assume  that  checking  if  an  ST  element  has  zero  ST  volume  is  equivalent  to 
checking  if  all  the  spatial  elements  associated  with  that  ST  element  have  zero 
volume.  Therefore,  for  this  purpose,  we  check  the  spatial-element  volumes. 

•  To  identify  the  spatial  elements  with  zero  volume,  which  should  have  zero  Jaco- 
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bian  at  all  the  integration  points,  instead  of  evaluating  the  Jacobians,  we  make 
the  determination  for  a  given  spatial  element  from  the  master-slave  relationship 
of  its  nodes. 


2.2  Design  of  the  master— slave  system 

The  data  we  need  to  provide  to  the  solver  is  quite  simple.  It  is  just  the  master-slave 
relationship  at  each  time  level.  However  there  are  some  restrictions,  and  here  we 
explain  the  three  that  we  want  to  emphasize. 

The  first  restriction  is  that  we  cannot  have  a  node  which  is  not  part  of  any  active 
(nonzero  volume)  spatial  element.  This  is  because  the  values  at  such  nodes  would  no 
longer  be  in  our  equation  system,  and  therefore  would  become  undefined.  If  because 
of  another  TC  such  a  node  comes  back  to  the  equation  system  later  as  part  of  an 
active  element,  it  would  add  an  undefined  component  to  the  equation  system. 

The  second  restriction  is  that  when  we  construct  the  ST  elements,  we  have  to  have 
matching  lateral  element-boundary  faces  between  the  active  adjacent  ST  elements. 
This  condition  cannot  be  checked  on  the  spatial  mesh.  Therefore  we  need  to  check  it 
on  the  ST  mesh. 

The  third  one  is  related  to  implementation.  The  master-slave  relationship  also 
extends  to  cases  when  we  have  boundary  conditions  on  the  master  and  slave  nodes. 
In  other  words,  the  conditions  at  the  master  node  also  apply  to  the  slave  nodes. 
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Chapter  3 

Conceptual  examples 


The  material  in  this  chapter  is  from  [59] . 


3.1  Contraction  and  expansion 


This  is  related  to  the  first  one  of  the  two  cases  of  TC  described  in  Chapter  2.  Con¬ 
traction  and  expansion  are  basically  the  same,  except  having  different  directions  in 
time  progression.  Figure  3.1  shows  a  contraction  example.  The  spatial  element  with 


Figure  3.1:  Contraction.  The  red  nodes,  3  and  5,  are  on  the  contraction  interface 
and  are  contacting.  The  white  nodes  are  the  slaves.  They  are  in  the  same  position 
as  their  masters,  but  for  visualization  purposes  we  slightly  shift  their  positions  in  the 
figure.  The  numbers  indicate  the  node  numbers  on  the  parent  mesh. 
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nodes  1  and  2,  for  example,  has  zero  volume  at  the  first  time  level.  However,  it  has 
nonzero  volume  at  the  second  time  level,  and  therefore  the  corresponding  ST  element 
has  nonzero  volume. 


3.2  Flapping 


This  is  related  to  the  second  one  of  the  two  cases  of  TC  described  in  Chapter  2. 


Figure  3.2  shows  the  red  and  blue  bars  at  three  instants  in  time  as  the  red  bar  crosses 


the  blue  bar.  Figure  3.3  shows,  for  the  flapping  motion,  the  ST  trajectories  of  the 
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Figure  3.2:  Flapping.  Red  and  blue  bars  at  different  instants  in  time  as  the  red  bar 
crosses  the  blue  bar. 


neighboring  ends  of  the  blue  and  red  bars.  Figure  3.4  shows  the  ST  element  edges 


1 - >  y 

Figure  3.3:  Flapping.  The  ST  trajectories  of  the  neighboring  ends  of  the  blue  and 
red  bars. 


for  the  line  separating  the  two  sides  of  the  domain  containing  the  blue  and  red  bars 
(shown  as  the  vertical  dashed  line  in  Figure  3.2). 
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Figure  3.4:  Flapping.  The  ST  element  edges  for  the  vertical  dashed  line  in  Figure  3.2. 

For  each  side  of  the  domain,  the  spatial  node  motions  along  the  ST  element  edges 
have  to  be  designed  in  a  fashion  that  does  not  lead  to  mesh  entanglement.  Figure  3.5 
shows  the  master-slave  relationship  for  the  blue-bar  side  of  the  domain,  and  Figure  3.6 
the  red-bar  side.  In  addition,  those  two  sides  are  in  a  master-slave  relationship  along 


Figure  3.5:  Flapping.  Blue-bar  side  of  the  ST  boundary  between  the  two  sides. 


the  vertical  dashed  line  in  Figure  3.2. 
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Figure  3.6:  Flapping.  Red-bar  side  of  the  ST  boundary  between  the  two  sides. 
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Chapter  4 


A  pair  of  symmetrically-flapping 
surfaces  (“Flapping  pair”) 


The  material  in  this  chapter  is  from  [59].  We  use  the  DSD/SST  method  in  the 
computations.  The  stabilization  parameter  Tgups  comes  from  the  tsupg  definition 
in  [67],  specifically  the  definition  given  by  Eqs.  (107)-(109)  in  [67],  which  can  also  be 
found  as  the  definition  given  by  Eqs.  (7)-(9)  in  [70],  with  /irgn  (=  ^rgnt)  and  r'LSic 
(=  Ulstc-hrgn)  given  by  Eqs.  (15)  and  (19)  in  [61].  In  solving  the  linear  equation 
systems  encountered  at  every  nonlinear  iteration  of  a  time  step,  the  GMRES  search 
technique  [42]  is  used  with  diagonal  preconditioner. 

4.1  Geometry,  motion  modeling  and  computational 
conditions 

The  pair  of  surfaces,  with  zero  thickness,  undergo  a  prescribed  sinusoidal  flapping,  as 
shown  in  Figures  4.1  and  4.2,  with  a  period  of  T  =  18.0  s.  The  projected  width  of 
the  flapping  surfaces  along  the  horizontal  axis  is  1.0  m.  The  density  and  kinematic 
viscosity  are  1.0  kg/m3  and  1.6667xl0-5  m2/s.  The  inflow  velocity  is  0.1  m/s.  The 
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Figure  4.1:  Flapping  pair.  Surface  positions  at  t  =  0  (also  18),  3,  6,  9,  12  and  15  s. 


Reynolds  number  based  on  these  length,  viscosity  and  velocity  scales  is  6,000. 

The  dimensions  of  the  computational  domain,  in  m,  are  40x20,  and  the  distance 
between  the  inflow  boundary  and  the  leading  edge  is  15  m.  The  boundary  conditions 
are  no-slip  on  the  flapping  surfaces,  uniform  horizontal  velocity  at  the  inflow,  zero- 
stress  at  the  outflow,  and  slip  at  the  upper  and  lower  boundaries.  We  tested  three 
different  meshes  to  see  the  influence  of  increased  refinement  in  space  and  time.  The 
meshes  have  a  structured  inner  zone  and  an  unstructured  outer  zone,  made  of  4- 
node  quadrilateral  and  3- node  triangular  elements,  respectively.  Table  4.1  shows,  for 
each  mesh,  the  number  nodes  along  each  flapping  surface,  number  of  nodes  in  the 
inner  zone,  total  number  of  nodes  and  elements,  and  the  number  of  time  steps  per 
flapping  cycle.  Figures  4. 3-4. 6  show  Mesh  1  and  the  inner  zones  for  all  three  meshes. 

During  the  flapping  motion,  only  the  inner  mesh  moves,  with  a  special,  algebraic 
mesh  moving  technique.  Figures  4. 7-4. 9  show,  for  Mesh  1,  the  inner  mesh  at  different 
instants  during  the  flapping  cycle. 
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Figure  4.2:  Flapping  pair.  Leading-edge  position,  velocity  and  acceleration  for  the 
upper  surface. 


Mesh 

Surface 

Nodes 

Inner 

Nodes 

Total 

Nodes 

Elements 

Time 

Steps 

i 

40 

28,379 

37,911 

47,628 

60 

2 

80 

32,439 

42,432 

52,590 

120 

3 

120 

36,499 

47,156 

57,958 

180 

Table  4.1:  Flapping  pair.  Mesh  data  and  the  number  of  time  steps  per  flapping 
cycle. 


4.2  Computations 


We  use  the  SUPS  version  of  the  DSD/SST  method  (see  [54,  56,  16]  for  the  termi¬ 
nology).  Prior  to  the  flapping  motion,  we  compute  2,500  time  steps  to  develop  the 
flow  field.  The  time  step  size  is  1.0  s,  with  3  nonlinear  iterations  per  time  step,  and 
the  corresponding  number  of  GMRES  iterations  are  50,  100  and  150.  The  flapping 
cycles  are  computed  with  5  nonlinear  iterations  per  time  step,  and  the  corresponding 
number  of  GMRES  iterations  are  150,  250,  450,  650  and  900.  We  note  that  since  it 
is  not  easy  to  show  that  the  solution  in  this  test  problem  is  accurate  in  terms  of  fluid 
physics,  we  decided  to  use  large  numbers  of  GMRES  iterations.  We  did  not  explore 
reducing  the  number  of  GMRES  iterations,  because  the  2D  computations  are  not 
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Figure  4.3:  Flapping  pair.  Mesh  1. 

costly.  We  compute  the  problem  for  two  complete  cycles  and  display  the  results  for 
the  second  flapping  cycle.  Figures  4.10-4.12  show,  for  the  three  meshes,  the  lift  and 
drag  experienced  by  the  flapping  surfaces. 

4.3  Results 

Figures  4.13-4.15  show,  for  Mesh  1,  the  velocity  magnitude  at  different  instants  during 
the  flapping  cycle. 
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Figure  4.4:  Flapping  pair.  Mesh  1.  Inner  zones 
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Figure  4.5:  Flapping  pair.  Mesh  2.  Inner  zones. 


35 


20 


Figure  4.6:  Flapping  pair.  Mesh  3.  Inner  zones. 
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Figure  4.7:  Flapping  pair.  Mesh  1  at  t  =  0  (also  18)  and  3  s. 
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Figure  4.8:  Flapping  pair.  Mesh  1  at  t  =  6  and  9  s. 
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Figure  4.9:  Flapping  pair.  Mesh  1  at  t  =  12  and  15  s. 
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—  Upper-Surface  Lift  —  Lower-Surface  Lift  Total  Drag 

Figure  4.10:  Flapping  pair.  Lift  and  drag  for  Mesh  1. 
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Figure  4.11:  Flapping  pair.  Lift  and  drag  for  Mesh  2. 
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—  Upper-Surface  Lift  —  Lower-Surface  Lift  Total  Drag 

Figure  4.12:  Flapping  pair.  Lift  and  for  Mesh  3. 
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0  0.1  0.2  0.3 

Figure  4.13:  Flapping  pair.  Velocity  magnitude  (in  m/s)  for  Mesh  1  at  t  =  0  and  3  s. 
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Figure  4.14:  Flapping  pair.  Velocity  magnitude  (in  m/s)  for  Mesh  1  at  t  =  6  and  9  s. 
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Figure  4.15:  Flapping  pair.  Velocity  magnitude  (in  m/s)  for  Mesh  1  at  t  =  12  and 
15  s. 
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Chapter  5 


Dragonfly  MAV 


The  material  in  this  chapter  is  from  [60] . 

5.1  Geometry  and  flapping-motion  modeling 

The  design  of  the  wings  is  similar  to  the  design  in  a  toy  MAV  [1],  The  body  is  the 
same  as  the  MAV  body  in  [50]. 

The  span  of  the  single  wing  is  46.7  mm  and  the  minimum,  maximum,  and  average 
chord  lengths  are  16.2,  19.2,  and  17.6  mm,  respectively  (see  Figure  5.1).  The  wings 
have  zero  thickness  and  undergo  prescribed  flapping,  as  shown  in  Figures  5.2  and  5.3, 
with  a  period  of  T  —  0.0365  s.  Figure  5.4  shows  the  position  of  the  leading-edge 
contact  point  over  time.  The  position  is  measured  from  the  body. 

The  density  and  kinematic  viscosity  are  1.225  kg/m3  and  1.461  xlO-5  m2/s.  The 
free-stream  velocity  is  4.5  m/s.  The  Reynolds  number  based  on  average  chord  length 
and  free-stream  velocity  is  5,423.  Three  cases  are  computed,  with  the  angle  of  attack 
a  =  0°,  5°,  and  10°. 

The  dimensions  of  the  computational  domain,  in  spans  of  a  single  wing,  are 
30  X  20  x20,  and  the  distance  between  the  inflow  boundary  and  the  leading  edge  is 
10  (see  Figure  5.5).  The  boundary  conditions  are  no-slip  on  the  wings  and  body 
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Figure  5.1:  Dragonfly  MAY.  Wing  dimensions. 


Figure  5.2:  Dragonfly  MAV.  Wing  configurations  at  t/T  =  0.0,  0.1,  0.2,  0.3,  0.4,  and 
0.5  (left  to  right  and  then  top  to  bottom). 
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Figure  5.3:  Dragonfly  MAY.  Wing  leading  edges  at  the  same  instants  as  in  Figure  5.2. 
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Figure  5.4:  Dragonfly  MAV.  Contact  point  position  along  the  leading  edge  over  a 
flapping  cycle. 
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of  the  MAV,  uniform  horizontal  velocity  at  the  inflow  boundary,  zero-stress  at  the 
outflow  boundary,  and  slip  at  the  upper,  lower,  and  side  boundaries. 

The  meshes  have  structured,  inner  zones  around  the  wings  and  an  unstructured, 
outer  zone.  Both  the  structured  and  unstructured  zones  consist  of  tetrahedral  ele¬ 
ments.  For  each  mesh,  Table  5.1  shows  the  number  of  nodes  and  elements.  Figure  5.6 
shows  top  view  of  the  wing  and  body  surface  meshes.  During  the  flapping  motion, 
only  the  mesh  in  the  inner  zones  move,  with  a  special,  algebraic  mesh  moving  tech¬ 
nique. 

The  structured,  inner  zones  consist  of  four  parts  corresponding  to  each  wing. 
Those  parts  each  have  3x2x2  structured  zones.  Each  zone  has  20x20x20  hexahedral 
clusters  made  of  6  tetrahedral  elements.  Figures  5.7  and  5.8  show,  for  a  =  0°,  the 
mesh  at  six  equally-spaced  instants  in  time  while  the  wings  are  closing.  The  zones 
between  the  upper  and  lower  wings  collapse  when  the  wings  close,  and  the  nodes  in 
the  neighboring  zones  also  collapse  accordingly.  We  note  that  a  wing  has  split  nodes 
except  on  the  leading  and  trailing  edges.  However,  when  the  wings  are  closed,  the 
nodes  on  the  upper  surface  of  the  upper  wings  and  the  lower  surface  of  the  lower 
wings  become  masters.  When  the  wings  are  partially  closed,  at  the  contact  point,  the 
nodes  on  the  lower  surface  of  the  upper  wing  are  also  masters  while  the  nodes  on  the 
upper  surface  of  the  lower  wings  are  slaves. 

5.2  Computational  conditions 

We  use  the  DSD/SST-SUPS  and  DSD/SST-VMST  (convective)  techniques  for  the 
first  two  and  last  two  nonlinear  iterations  of  each  time  step.  The  stabilization  pa¬ 
rameter  tsups  comes  from  the  tsupg  definition  in  [67],  specifically  the  definition  given 
by  Eqs.  (107)-(109)  in  [67],  which  can  also  be  found  as  the  definition  given  by  Eqs. 
(7)-(9)  in  [70],  with  z/Lsic  from  Eq.  (17)  in  [70].  The  time-step  size  is  4.51xl0-4  s. 
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Figure  5.5:  Dragonfly  MAV.  Computational  domain  and  mesh  setup.  Outer  bound¬ 
aries  (gray),  boundaries  of  the  inner,  structured  meshes  (blue)  and  body  (green). 


Figure  5.6:  Dragonfly  MAV.  Surface  mesh  at  t/T  =  0.5. 
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Figure  5.7:  Dragonfly  MAV.  Mesh  (cut  mid-chord)  at  the  same  instants  as  in  Fig¬ 
ure  5.2  (left  to  right  then  top  to  bottom). 
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Figure  5.8:  Dragonfly  MAV.  Mesh  (cut  mid-span)  at  the  same  instants  as  in  Figure  5.2 
(left  to  right  then  top  to  bottom). 
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Table  5.1:  Dragonfly  MAV.  Number  of  nodes  (nn)  and  elements  (ne)  in  the  meshes 
used. 


Interface  Mesh 

Single  Wing 

nn 

802 

ne 

1,600 

Body 

nn 

8,521 

ne 

16,630 

Volume  Mesh 

Inner 

nn 

405,002 

ne 

2,303,920 

Full  (a  =  0°) 

nn 

1,143,613 

ne 

6,844,706 

Full  (a  =  5°) 

nn 

1,227,618 

ne 

7,353,540 

Full  (a  =  10°) 

nn 

1,152,367 

ne 

6,896,762 

For  each  angle  of  attack,  prior  to  the  flapping  motion,  we  compute  550  time  steps 
with  the  geometry  at  t  =  0  to  develop  the  flow  field.  For  the  first  500  time  steps,  only 
half  of  the  computational  domain  is  used  and  a  slip  boundary  condition  is  enforced 
on  the  symmetry  plane.  The  results  are  then  copied  to  the  other  half  of  the  mesh  for 
the  final  50  time  steps  of  flow  field  development.  The  inflow  velocity  of  4.5  m/s  is 
reached  by  a  sinusoidal  ramping  over  the  first  150  time  steps,  starting  from  0.0  m/s. 
In  computing  the  developed  flow  field,  the  number  of  GMRES  iterations  per  nonlinear 
iteration  is  150,  350,  450,  and  800.  In  computing  the  flapping  cycles,  the  number  of 
GMRES  iterations  is  250,  500,  750,  and  1,000.  We  compute  three  flapping  cycles  and 
display  the  results  for  the  third  cycle. 


5.3  Results 


We  first  present  (in  Figures  5.9-5.14),  only  for  a  =  10°,  results  over  (or  in  relationship 
to)  the  MAV  body  and  wing  surfaces.  Figure  5.9  shows  the  helicity  isosurfaces.  The 
flow  field  near  the  wings  is  almost  symmetric,  but  the  flow  behind  the  MAV  is  not. 
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Figures  5.10  and  5.11  show  pressure  on  the  body  and  the  wing  surfaces.  The  pressure 
is  almost  symmetric,  and  therefore  we  use  the  left  and  right  sides  of  the  wing  pictures 
for  the  upper  and  lower  wings.  For  the  body,  however,  both  sides  show  the  upper 
surface.  Figures  5.12  and  5.13  show  the  magnitude  of  shear  stress  on  the  body  and 
wing  surfaces.  Again,  we  use  the  left  and  right  sides  of  the  wing  pictures  for  the  upper 
and  lower  wings,  and  both  sides  of  the  body  for  the  upper  surface.  Figure  5.14  shows 
the  pressure  difference  between  the  upper  and  lower  surfaces.  The  left  side  is  for  the 
upper  wing,  and  the  right  side  for  the  lower  wing.  For  the  closed  parts  of  the  wings, 
both  sides  show  the  difference  between  the  lower  surface  of  the  lower  wing  and  the 
upper  surface  of  the  upper  wing.  Lift  and  drag  forces  are  shown  in  Figures  5.15-5.17. 
The  forces  are  separated  into  the  upper  and  lower  wings,  the  closed  wings,  and  the 
body. 
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Figure  5.9:  Dragonfly  MAV.  Helicity  isosurfaces  (±5  and  ±10  m2/s2)  for  a  =  10°  at 
t/T  =  0.0,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  and  0.9  (left  to  right  and  then  top  to 
bottom).  Blue  is  for  negative  values,  and  red  for  positive. 
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Figure  5.10:  Dragonfly  MAV.  Pressure  (Pa)  for  a  =  10°  on  the  body  and  wings  at 
the  same  instants  as  in  Figure  5.9.  The  upper  surface  of  the  upper  wing  (left  side) 
and  the  lower  surface  of  the  lower  wing  (right  side). 


55 


40 


-100  -50  0  50  100 


Figure  5.11:  Dragonfly  MAV.  Pressure  (Pa)  for  a  =  10°  on  the  body  and  wings  at 
the  same  instants  as  in  Figure  5.9.  The  lower  surface  of  the  upper  wing  (left  side) 
and  the  upper  surface  of  the  lower  wing  (right  side) .  The  white  regions  are  the  closed 
parts  of  the  wings. 
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Figure  5.12:  Dragonfly  MAV.  Magnitude  of  the  shear  stress  (Pa)  on  the  body  and  the 
wing  surfaces  at  the  same  instants  as  in  Figure  5.9.  The  upper  surface  of  the  upper 
wing  (left  side)  and  the  lower  surface  of  the  lower  wing  (right  side). 
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Figure  5.13:  Dragonfly  MAV.  Magnitude  of  the  shear  stress  (Pa)  on  the  body  and 
the  wing  surfaces  at  the  same  instants  as  in  Figure  5.9.  The  lower  surface  of  the 
upper  wing  (left  side)  and  the  upper  surface  of  the  lower  wing  (right  side) .  The  white 
regions  are  the  closed  parts  of  the  wings. 
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Figure  5.14:  Dragonfly  MAV.  Pressure  difference  (Pa)  between  the  lower  and  upper 
surfaces  for  a  =  10°  at  the  same  instants  as  in  Figure  5.9.  For  the  upper  wing  and 
closed  wings  (left  side)  and  for  the  lower  wing  and  closed  wings  (right  side). 
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—  Upper  Wings  —  Lower  Wings  —  Closed  Wings 


Body  — Total 
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—  Upper  Wings  —  Lower  Wings  —  Closed  Wings 
Body  — Total 

Figure  5.15:  Dragonfly  MAY.  Lift  (top)  and  drag  (bottom)  for  a  =  0°. 
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—  Upper  Wings  —  Lower  Wings  —  Closed  Wings 
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—  Upper  Wings  —  Lower  Wings  —  Closed  Wings 
Body  — Total 

Figure  5.16:  Dragonfly  MAY.  Lift  (top)  and  drag  (bottom)  for  a  =  5°. 
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—  Upper  Wings  —  Lower  Wings  —  Closed  Wings 


Body  — Total 
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—  Upper  Wings  —  Lower  Wings  —  Closed  Wings 
Body  — Total 

Figure  5.17:  Dragonfly  MAY.  Lift  (top)  and  drag  (bottom)  for  a  =  10°. 
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Chapter  6 

Aortic  valve  with  coronary  arteries 

The  material  in  this  chapter  is  from  [58]. 

6.1  Geometry  and  motion  modeling 

We  create  a  typical  aortic- valve  model  based  on  pictures,  such  as  the  one  in  [21].  The 
model,  shown  in  Figures  6.1  and  6.2,  has  three  leaflets  with  two  outlets,  corresponding 
to  coronary  arteries,  and  one  main  outlet,  corresponding  to  the  beginning  of  the  aorta. 
The  outlets  are  extended  straight  in  each  direction.  The  bulges  are  called  sinuses. 
The  left  and  right  coronary  arteries  are  attached  to  the  left  and  right  aortic  sinuses, 
respectively,  and  the  other  sinus  is  called  the  posterior  aortic  sinus. 

The  inlet  and  main  outlet  diameters  are  23  mm,  which  correspond  to  a  typical 
aorta.  The  two  coronary  artery  diameters  are  2.9  mm. 

We  prescribe  the  motion  of  the  three  leaflets  so  that  the  valve  closes  and  opens 
with  a  period  of  T  =  0.6  s.  The  inflow  velocity  is  specified  such  that  the  average 
flow  rate  is  5,000  ml/min.  The  flow  rate  is  a  time- variant  function  of  the  horizontal 
projection  of  the  open  mouth  area  of  the  valve  when  it  is  open,  and  of  the  inlet-side 
volume  change  while  the  valve  is  changing  shape  after  it  closes  (see  Figure  6.3).  The 
density  and  kinematic  viscosity  of  the  blood  are  1,000  kg/m3  and  4.0xl0-6  m2/s. 
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Figure  6.1:  Aortic  valve  with  coronary  arteries.  Model  geometry.  Aorta,  leaflets, 
sinuses,  and  coronary  arteries.  The  left  coronary  artery  is  on  the  right  in  the  figure, 
and  the  right  coronary  artery  is  on  the  left  in  the  figure. 
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Figure  6.2:  Aortic  valve  with  coronary  arteries.  Leaflets  at  t/T  =  0.0,  0.1,  0.2,  0.3, 
0.4,  0.5,  0.6,  0.7,  0.8,  and  0.9  (left  to  right  and  then  top  to  bottom). 
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Figure  6.3:  Aortic  valve  with  coronary  arteries.  Inflow  profile. 


The  boundary  conditions  are  no-slip  on  the  arterial  walls  and  valves,  traction- free 
at  the  outflow  boundaries,  and  uniform  velocity  at  the  inflow  boundary.  The  surface 
mesh  on  the  leaflets  (shown  in  Figure  6.4)  is  made  up  of  8,448  nodes  and  16,440 
triangular  elements.  The  arterial- wall  surface  mesh  is  shown  in  Figure  6.5.  The  mesh 
has  structured,  inner  zones  around  the  leaflets  and  an  unstructured,  outer  zone.  The 
inner  zones  consist  of  tetrahedral,  pyramid-shaped,  and  wedge-shaped  elements,  and 
the  outer  zone  consists  of  tetrahedral  elements.  The  volume  mesh  is  made  up  of 
1,417,910  nodes  and  4,184,614  mixed  elements.  The  mesh  near  the  valve  is  shown  in 
Figure  6.6. 

During  the  prescribed  motion,  only  the  inner  zones  move  with  a  special,  algebraic 
mesh-moving  technique.  The  positions  of  the  nodes  in  the  inner  zones  are  created  by 
linearly  interpolating  the  surface  mesh  of  the  zones  from  the  closed  position  to  the 
open  position.  Seventy-nine  layers  of  nodes  are  extruded  from  both  the  upper  and 
lower  surfaces  of  the  leaflets. 

When  the  valve  is  completely  open,  all  of  the  fluid  nodes  extruded  from  the  upper 
surface  are  slaves  to  the  upper  surface,  and  all  of  the  nodes  extruded  from  the  lower 
surface  are  masters.  As  the  valve  closes,  it  leaves  one  layer  of  nodes  attached  to  the 
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upper  surface  behind,  and  the  leaflet  surface  moves  to  the  next  node  position,  making 
those  nodes  slaves  to  the  lower  surface.  When  the  valve  is  completely  closed,  all  of 
the  fluid  nodes  extruded  from  the  upper  surface  are  masters  and  all  of  the  nodes 
extruded  from  the  lower  surface  are  slaves  to  the  lower  surface. 


6.2  Computational  conditions 

We  use  the  DSD/SST-SUPS  and  DSD/SST-VMST  (convective)  techniques  for  the 
first  two  and  last  two  nonlinear  iterations  of  each  time  step.  The  stabilization  pa¬ 
rameter  rgups  comes  from  the  tsupg  definition  in  [67] ,  specifically  the  definition  given 
by  Eqs.  (107)-(109)  in  [67],  which  can  also  be  found  as  the  definition  given  by  Eqs. 
(7)-(9)  in  [70],  with  nLsic  from  Eq.  (17)  in  [70].  Prior  to  the  prescribed  motion, 
we  compute  150  time  steps  with  the  geometry  at  t  =  0  to  develop  the  flow  field. 
The  viscosity  of  4.0  xlO-6  m2/s  is  reached  by  ramping  over  the  first  50  time  steps 
starting  from  the  viscosity  1.31  xlO-3  m2/s.  The  ramping  profile  for  the  viscosity  is 
designed  to  result  in  a  linear  ramping  for  the  Reynolds  number.  The  time-step  size  is 
6.33xl0-3  s  during  flow-field  development,  and  3.00xl0-3  s  for  the  prescribed-motion 
cycles.  In  computing  the  developed-flow  field,  the  number  of  GMRES  iterations  per 
nonlinear  iteration  is  150,  350,  450,  and  800.  In  computing  the  flapping  cycles,  the 
number  of  GMRES  iterations  is  250,  500,  750,  and  1,000.  We  compute  two  cycles 
and  display  the  results  for  the  second  cycle. 

6.3  Results 

The  global  mass-balance  error,  normalized  by  the  average  flow  rate,  is  less  than 
10%.  Figure  6.7  shows  the  flow  rate  at  the  outlets  of  the  coronary  arteries.  There 
is  some  “negative  outflow”  (i.e.  inflow)  from  the  coronary  arteries,  however,  the  wall 
shear  stress  (WSS)  on  the  long  pipes  are  large  enough  to  stabilize  the  overall  system. 
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Figure  6.4:  Aortic  valve  with  coronary  arteries.  The  leaflets  surface  mesh  at  the  same 
instants  as  in  Figure  6.2. 
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Figure  6.5:  Aortic  valve  with  coronary  arteries.  Mesh  of  the  aortic  valve,  sinuses,  and 
coronary  arteries. 
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Figure  6.6:  Aortic  valve  with  coronary  arteries.  Mesh  around  the  leaflets  at  t/T 
0.0,  0.2,  0.4,  and  0.6. 
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t/T 

Figure  6.7:  Aortic  valve  with  coronary  arteries.  Flow  rate  at  the  outlets  of  the 
coronary  arteries. 

Figures  6.8-6.10  show  a  volume  rendering  of  the  velocity  magnitude.  Nonsymmetric 
and  complex  flow  patterns  are  observed  behind  the  valve.  Figures  6.11  and  6.12 
show  the  velocity  magnitude  on  the  “coronary  plane,”  where  the  coronary  arteries 
are  connected  to  the  aorta,  and  the  “above-sinus  plane,”  which  is  18  mm  downstream 
from  the  coronary  plane.  The  mainstream  flow  oscillates  away  from  the  sinuses.  This 
is  mainly  due  to  the  jet  from  the  contact  of  the  leaflets  of  the  sinuses. 

Figures  6.13  and  6.14  show  pressure  difference  between  the  lower  and  upper  sur¬ 
faces  of  the  leaflets.  We  exclude  the  parts  where  the  leaflets  are  in  contact.  Fig¬ 
ures  6.15-6.18  show  WSS  on  the  leaflet  surfaces.  The  WSS  on  the  lower  surfaces  of 
the  three  leaflets  are  somewhat  similar  to  each  other.  However,  on  the  upper  surfaces, 
the  WSS  for  the  leaflets  of  the  coronary  sinuses  are  different  from  the  WSS  for  the 
leaflet  of  the  posterior  sinus.  Figure  6.19  shows  OSI  on  the  leaflet  surfaces.  The  WSS 
vector  is  projected  onto  the  open  configuration  to  calculate  OSI  (see  [52]). 
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Figure  6.8:  Aortic  valve  with  coronary  arteries.  Volume  rendering  of  the  velocity 
magnitude  (m/s)  at  t/T  =  0.0,  0.1,  0.2,  and  0.3. 
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Figure  6.9:  Aortic  valve  with  coronary  arteries.  Volume  rendering  of  the  velocity 
magnitude  (m/s)  at  t/T  =  0.4,  0.5,  0.6,  and  0.7. 
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Figure  6.10:  Aortic  valve  with  coronary  arteries.  Volume  rendering  of  the  velocity 
magnitude  (m/s)  at  t/T  =  0.8,  and  0.9. 
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Figure  6.11:  Aortic  valve  with  coronary  arteries.  The  velocity  magnitude  on  the 
coronary  plane  (m/s)  at  the  same  time  instants  as  in  Figure  6.2. 
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Figure  6.12:  Aortic  valve  with  coronary  arteries.  The  velocity  magnitude  (m/s)  on 
the  above- sinus  plane  at  the  same  time  instants  as  in  Figure  6.2. 
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Figure  6.13:  Aortic  valve  with  coronary  arteries.  Pressure  difference  (Pa)  between 
the  lower  and  upper  surfaces  of  the  three  leaflets  at  t/T  =  0.0,  0.1,  0.2,  0.3,  0.4,  and 
0.5. 
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Figure  6.14:  Aortic  valve  with  coronary  arteries.  Pressure  difference  (Pa)  between 
the  lower  and  upper  surfaces  of  the  three  leaflets  at  t/T  =  0.6,  0.8,  0.8,  and  0.9. 
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Figure  6.15:  Aortic  valve  with  coronary  arteries.  WSS  (Pa)  on  the  lower  surface  of 
the  three  leaflets  at  the  same  time  instants  as  in  Figure  6.13. 
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Figure  6.16:  Aortic  valve  with  coronary  arteries.  WSS  (Pa)  on  the  lower  surface  of 
the  three  leaflets  at  the  same  time  instants  as  in  Figure  6.14. 
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Figure  6.17:  Aortic  valve  with  coronary  arteries.  WSS  (Pa)  on  the  upper  surface  of 
the  three  leaflets  at  the  same  time  instants  as  in  Figure  6.13. 
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Figure  6.18:  Aortic  valve  with  coronary  arteries.  WSS  (Pa)  on  the  upper  surface  of 
the  three  leaflets  at  the  same  time  instants  as  in  Figure  6.14. 
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Figure  6.19:  Aortic  valve  with  coronary  arteries.  OSI  on  the  lower  (left)  and  upper 
(right)  surfaces  of  the  leaflets.  The  leaflets  are  in  the  fully  open  configuration,  and 
left,  right,  and  posterior  aortic  sinuses  from  top  to  bottom,  respectively. 
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Chapter  7 

Concluding  Remarks 


The  material  in  this  chapter  is  from  [59].  We  have  presented  an  interface-tracking 
(moving-mesh)  method  that  addresses  the  computational  challenges  associated  with 
contact  between  moving  solid  surfaces  and  other  cases  of  TC,  including  those  in  car¬ 
diovascular  FSI,  parachute  FSI,  and  flapping- wing  aerodynamics,  with  the  primary 
target  in  cardiovascular  FSI  being  heart  valves.  It  is  a  new  version  of  the  DSD/SST 
method,  and  we  call  it  ST-TC.  It  is  based  on  special  mesh  generation  and  update,  and 
a  master-slave  system  that  maintains  the  connectivity  of  the  parent  mesh  when  there 
is  contact  or  other  TC.  This  makes  the  method  an  efficient,  practical  alternative  to 
using  unstructured  ST  meshes,  but  without  giving  up  on  the  accurate  representation 
of  the  interface  or  consistent  representation  of  the  interface  motion,  even  where  there 
is  contact  or  other  TC.  We  explained  the  method  with  conceptual  examples,  and  pre¬ 
sented  successful  2D  and  3D  computations  with  models  representative  of  the  classes 
of  problems  we  are  targeting.  We  are  comfortable  with  concluding  that  the  ST-TC 
method  has  the  interface-tracking  accuracy,  the  TC  flexibility,  and  the  computational 
practicality. 
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